\documentclass[10pt,a4paper]{article} 

\usepackage{ctex} % 中文支持
\usepackage[top=2.5cm, bottom=2.5cm, left=2.5cm, right=2.5cm]{geometry} % 页边距
\usepackage{amsmath, amssymb} % 数学公式与符号
\usepackage{graphicx}
\usepackage{pythonhighlight}
\usepackage{url} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{titling}
\setlength{\droptitle}{-2cm} % 标题上移

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%文档的题目、作者与日期
\author{五六七 }
\title{洛伦兹模型与混沌现象 }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\maketitle

\begin{abstract}
研究洛伦兹方程与混沌现象。
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setcounter{tocdepth}{2}
%\renewcommand\contentsname{目录}
%
%\renewcommand {\baselinestretch} {1.3}\normalsize 
%\tableofcontents 
%\renewcommand {\baselinestretch} {1.0}\normalsize


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{问题描述}

美国气象学家洛伦兹在研究大气运动时，简化大气对流模型，得到如下一阶自治乘微分方程组，
\begin{eqnarray}
\left\{\begin{array}{rcl}
\frac{dx}{dt} &=& \sigma(y-x), \\ 
\frac{dy}{dt} &=& \rho x - y -xz, \\ 
\frac{dz}{dt} &=& xy - \beta z, \\ 
\end{array}\right. 
\end{eqnarray}
其中 $\sigma,\rho,\beta$ 是三个参数。
研究这三个参数的不同选择，判断系统是否会产生混沌现象。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section{建立模型}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{编程计算A}
首先载入数值计算的numpy模块，科学计算的scipy模块里的求解微分方程的函数odeint, 以及画图的pylab模块。
\begin{python}
from scipy.integrate import odeint
import numpy as np
import pylab as plt
\end{python}

为了进行一致性比较，每次运行取相同随机数，为此设置随机数种子。
设置洛伦兹模型的三个参数的值。
最后使用lambda函数定义微分方程组的右端项，一共有三个未知函数。
\begin{python}
np.random.seed(2)  
sigma=10; rho=28; beta=8/3;
g=lambda f,t: [sigma*(f[1]-f[0]), rho*f[0]-f[1]-f[0]*f[2], f[0]*f[1]-beta*f[2]]   
\end{python}

用随机数设置初始值，定义自变量（时间变量）的取值区间和分点，然后调用odeint函数求洛伦兹方程的数值解。
\begin{python}
s01=np.random.rand(3)
t0=np.linspace(0,40,2000)
s1=odeint(g,s01,t0)
\end{python}

改变初始值一点点，再求微分方程的数值解。
\begin{python}
s02=s01+0.000001
s2 = odeint(g,s02,t0)  
\end{python}

书里的这两行代码可以不用。
\begin{python}
#plt.rc('text',usetex=True)
#plt.rc('font',size=16)
\end{python}

画出第一个初始值的轨线。
\begin{python}
ax=plt.subplot(111, projection='3d')
plt.plot(s1[:,0],s1[:,1],s1[:,2],'b-')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
\end{python}

\begin{figure}[ht!]
\centering
\includegraphics[height=8cm, width=12cm]{lorenz01.png}
\caption{第一个初始值的洛伦兹曲线 }
\end{figure}

画出第二个初始值的轨线。
\begin{python}
bx=plt.subplot(111, projection='3d')
bx.plot(s2[:,0],s2[:,1],s2[:,2],'b-')
bx.set_xlabel('$x$')
bx.set_ylabel('$y$')
bx.set_zlabel('$z$')
\end{python}

\begin{figure}[ht!]
\centering
\includegraphics[height=8cm, width=12cm]{lorenz02.png}
\caption{第二个初始值的洛伦兹曲线 }
\end{figure}

画出两个初始值的轨线的差，如图3所示。
\begin{python}
cx=plt.subplot(111)
cx.plot(t0,s1[:,0]-s2[:,0],'-')
cx.set_xlabel('$t$')
cx.set_ylabel('$x_1(t)-x_2(t)$',rotation=90)
\end{python}

\begin{figure}[ht!]
\centering
\includegraphics[height=8cm, width=12cm]{lorenz03.png}
\caption{微小差异的初始值导致的洛伦兹曲线的差异 }
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{编程计算B}

使用scipy.integrate的odeint函数求解微分方程。
使用 mplot3d 子模块的 Axes3D 对象来创建三维空间里的曲线。
注意 mplot3d是以matplotlib的基础。
\begin{python}
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
\end{python}

定义洛伦兹微分方程组的右端的三个函数。
\begin{python}
def lorenz(w,t):
#    sigma=10; rho=28; beta=8/3
    sigma=5; rho=20; beta=3
    x,y,z=w
    return np.array([sigma*(y-x),rho*x-y-x*z,x*y-beta*z])
\end{python}

使用 scipy 模块的 odeint 函数来求解常微分方程组，第一个参数是洛伦兹方程，第二个参数是初始位置，第三个参数是时间的离散值。求解的结果保存在变量 solution 里。
\begin{python}
t=np.arange(0,30,0.01)
solution=odeint(lorenz,[0,1,0],t)
\end{python}

画出轨线图。变量 solution里存放的三个分量分别是轨线的x,y,z坐标。
\begin{python}
fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')
ax.plot(solution[:,0],solution[:,1],solution[:,2],'r--')
plt.savefig('lorenz04.png')
\end{python}

\begin{figure}[ht!]
\centering
\includegraphics[height=8cm, width=12cm]{lorenz04.png}
\caption{没有出现混沌现象的洛伦兹曲线 }
\end{figure}

图像如图4所示。注意到三个参数的取值， $\sigma=5, \rho=20, \beta=3$ 的图像，与 $\sigma=10, \rho=28, \beta=8/3$ 的图像有较大差别。最后可使用 type 函数来查看变量 fig 和变量 ax 分别是什么对象。
\begin{python}
type(fig)    
type(ax)
\end{python}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{回答问题}

我们选取洛伦兹模型的三个参数的取值，以及两个有着细微差别的初始值，观察到轨线在某个时刻之后产生较大差别，从而出验证了混沌现象。然后选取三个参数的另一组取值，没有观察到轨线有不同的中心。


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{参考文献 }
\begin{thebibliography}{99}

%\bibitem{dingtongren} 丁同仁、李承治，常微分方程教程，高等教育出版社，2022年3月第三版。
\bibitem{sishoukui-2} 司守奎,孙玺菁. \emph{Python数学建模算法与应用}, 国防工业出版社. 2022年1月第1版. 
\bibitem{lorenz} Lorenz, E. N. \emph{Deterministic Nonperiodic Flow}. Journal of the Atmospheric Sciences, 20, 130-141, 1963.

\bibitem{wolfram} Wolfram MathWorld.  \url{https://mathworld.wolfram.com/LorenzAttractor.html}

\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





